/**
 * Everything coarse graining, only good for d=2.
 */

public class Grain
{
	Ising ic;
	protected int L_reduction;
	private int grainL;
	private int grainN;
	private int [] g; // coarse grained lattice
	private final int forwardSteps = 2; // # of steps before coarse graining

	public Grain (Ising ic, int L_reduction)
	{
		this.ic = ic;
		this.L_reduction = L_reduction;
		if (ic.L%L_reduction != 0) 
		{
			System.err.println ("L does not divide coarse grain L.");
			System.exit (1);
		}
		grainL = ic.L / L_reduction;
		grainN = ic.N; for (int i=0; i<ic.d; ++i) grainN /= L_reduction;
		g = new int [grainN];
	}

	public Grain (Ising ic) { this (ic, 4); }

	private void coarseGrain ()
	{
		for (int i=0; i<grainN; ++i) g[i] = 0; // reset
		for (int i=0; i<ic.N; ++i)
		{
            int L_mult = 1;
            int G_mult = 1;
            int index = 0;
            for (int j=0; j<ic.d; ++j)
            {
                index += G_mult*(((i/L_mult)%ic.L)/L_reduction);
                L_mult *= ic.L;
                G_mult *= grainL;
            }
			g[index] += ic.getSpin (i);
		}
	}

	/**
	 * Estimate center by finding center of biggest grain.
	 */
	public int estCenter () 
	{
		coarseGrain ();
		int minGrain = 0;
		int [] x0 = new int [ic.d];
		int L_mult = 1;
		int center = 0;

		boolean [] s = ic.getLattice (); // back up lattice
		for (int i=0; i<forwardSteps; ++i) ic.update ();

		for (int i=0; i<grainN; ++i) if (g[i] < g[minGrain]) minGrain = i;
		for (int i=0; i<ic.d; ++i)
		{
			x0[i] = (minGrain % grainL)*L_reduction + L_reduction/2; 
			minGrain /= grainL;
		}
		for (int i=0; i<ic.d; ++i)
		{
			center += x0[i]*L_mult;
			L_mult *= ic.L;
		}
		// restore lattice
		for (int i=0; i<ic.N; ++i) { ic.setSpin (i, s[i] ? 1 : -1); }
		return center;
	}

	/**
	 * Distance of target from reference accounting for periodic boundaries.
	 * @param ref the reference location on one axis
	 * @param target the target location on one axis
	 */
	protected int delta (int ref, int target)
	{
		int dist;
		if (Math.abs(target - ref) > ic.L/2) // wrap
			dist = (target<ref) ? (target+ic.L-ref) : -(ref+ic.L-target);
		else // no wrap
			dist = target - ref;
		return dist;
	}

	protected double realDistance (int center, int site)
	{
		int [] x0 = new int [ic.d]; // distance from center
		int L_mult = 1;
		for (int j=0; j<ic.d; ++j)
		{
			// absolute value redundant in magnitude calculation
			x0[j] = Math.abs(delta ((center/L_mult)%ic.L, (site/L_mult)%ic.L));
			L_mult *= ic.L;
		}
		return Utils.magnitude (x0);
	}

	public int findCenter ()
	{
		int ref = estCenter ();
		int [] x0 = new int [ic.d]; // Average distance from reference
		int num = 0;
		int L_mult = 1;
		int center = 0;
		for (int i=0; i<ic.N; ++i)
		{
		    /*	
            // check to make sure the target spin is within a certain 
			// hypercube radius of the center estimate (only check adjacent 
			// coarse grain blocks
			boolean consider = true;
			L_mult = 1;
			for (int j=0; j<ic.d; ++j)
			{
				if (((ref/L_mult-i/L_mult)%ic.L+ic.L)%ic.L > 1.5*L_reduction)
					consider = false;
				L_mult *= ic.L;
			}
			if (!consider) continue; // skip if we don't consider
            */
            if (realDistance (ref, i) > ic.L/4) continue;

			if (ic.getSpin (i) == -1) // Lined up with field, always down
			{
				++num; // For center of mass calculation
				L_mult = 1;
				for (int j=0; j<ic.d; ++j)
				{
					x0[j] += delta ((ref/L_mult)%ic.L, (i/L_mult)%ic.L);
					L_mult *= ic.L;
				}
			}
		}
		for (int i=0; i<ic.d; ++i) 
			if (num > 0) // to prevent divide by zero
				x0[i] /= num; // normalize

		L_mult = 1;
		for (int i=0; i<ic.d; ++i)
		{
			center += ((ref/L_mult)%ic.L + x0[i])*L_mult;
			L_mult *= ic.L;
		}
		return (center + ic.N) % ic.N;
	}
}
